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Electrons in double-layer quantum well systems behave like pseudo-spin 1/2 particles 
where the up and down "spin" represent localized states in each of the layers. The magnet- 
ically induced Wigner crystals in these systems are therefore crystals of these pseudo-spin 
1/2 particles. We have calculated the phase diagram of the bilayer Wigner crystals using 

> 

a variational scheme which explores a continuum of lattice and spin structure. Five stable 

O 

CO 
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l/-) ' encounters the following sequence of transitions as the filling factor is increased from zero, 
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(the same sequence also occurs if one increases the "effective" layer separation starting from 

d ; 

S zero, with tunneling strength and filling factor held fixed) : 

O 
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crystal phases are found. For given tunneling strength and layer separation, one typically 



(I) (One-component) hexagonal structure — > (II) Centered rectangular structure 
(III) Centered square structure — > (IV) Centered rhombic structure 



— ► (V) Staggered hexagonal structure. 



Crystal I is a ferromagnet in pseudo-spin space. All other crystals (II to V) have mixed 
ferro- and antiferromagnetic orders, which are generated by layer tunneling and interlayer 
repulsion respectively. The relative strength of these two magnetic orders vary continuously 
with external parameters, (i.e. the ratio of layer separation to magnetic length, the tunneling 
gap to Coulomb interaction, etc). The lattice structures I, III, V are "rigid" whereas II and 
V are "soft" , in the sense that the latter two vary with external parameters and the former 
three do not. Another important feature of the phase diagram is the existence a multicritical 
point and a critical end-point, which allows all crystals (except V) to transform into one 
another continuously. While our findings are based on a variational calculation, one can 



conclude on physical grounds that the mixed ferromagnetic-antiferromagnetic order as well 
as the pseudospin-lattice coupling should be general features of most bilayer Wigner crystals. 

PACS no. 75.10.-b, 73.20.Dx, 64.60.Cn 
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I. INTRODUCTION 



Recent experiments |]J have shown that bilayer quantum well systems in strong magnetic 
fields behave like quantum Hall systems with (pseudo-)spin 1/2, where the up and down 
pseudo-spins correspond to electron states residing in one or the other of the layers. The 
novelty of this new quantum Hall system is demonstrated by the quantum Hall effect at 
filling factors v = 1/2 and v — 1, which corresponds to each layer being 1/4 filled and 1/2 
filled respectively ||T| . Since even denominator quantum Hall fluids do not exist in fully spin 
polarized single-layer systems, these states can only result from the correlations between 
electrons in different layers. 

One major difference between single-layer and bilayer systems is their energy scales. 
In the single-layer case, if one restricts to the lowest Landau level, the only remaining 
energy scale is the Coulomb interaction between electrons e 2 \fn^ where n is the electron 
density, (n is related to the filling factor v as n = z//(27r£ 2 ), where £ = (ch/eB) 1 ^ 2 is the 
magnetic length.) For bilayer systems, there are three energy scales: (a) the tunneling 
energy A between the layers, (precise definition given later), (b) the inter-layer Coulomb 
interaction e 2 /D, where D is the separation between the layers and, (c) the intra-layer 
Coulomb interaction e 2 \/n, where n = u/(2tt£ 2 ) is the total electron density of both layers. 
Because of these energy scales, the system can fall into different physical regimes depending 
on their relative ratios. The important physical regimes are: 

(i) Two single-layer regime, e 2 ^/n ^> e 2 /D,A ; 

(ii) Two-component (or correlation) regime, e 2 y/n ~ e 2 /D 3> A ; 

(iii) One-component regime, e 2 ^/n ~ e 2 /D < A. 

Regime (i) corresponds to the case of large layer spacing. In this case, the systems reduces 
to two weakly interacting single-layer systems. Regime (iii) corresponds to sufficiently small 
layer spacing so that the effects of tunneling becomes important. In this case, the electrons 
lie in the "symmetric" state of the quantum well. The system reduces again to the single- 
layer case. Regime (ii) is very different. Despite the weakness of layer tunneling, electrons 
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are strongly coupled to each other through the Coulomb interaction. The newly discovered 
v = 1/2 and v — 1 quantum Hall states are found in this regime In particular, the v = 
1/m states (m odd) are believed to behave like quantum ferromagnets 0. [ By one and two- 
component, we have in mind a pseudo-spin 1/2 representation of the electron wavefunction. 
We shall discuss this representation shortly. ] 

For single-layer systems, it is well known that the quantum Hall fluids eventually lose 
their stability at sufficiently low filling factors to the Wigner crystals, which are essentially 
products of electron Gaussians arranged on a lattice, properly antisymmetrized to satisfy 
the Pauli principle. The emergence of Wigner crystals at low filling is inevitable, as electron 
exchange has weakened so much, that the correlation energy of the fluid falls below the 
Madelung energy of the Wigner solid. In the last few years, many experiments (transport 
H 0, threshold field ||, magnetophonon 0, and luminlescence ||) have indicated the 
existence of crystalline characteristics at various filling factors, some are as high as around 
v = 1/5 J4j] ||. While there are questions as to whether the insulating states around 
v = 1/5 are Wigner crystals or other kinds of novel insulators [[J, the general view is that 
the insulating states at lower filling factors will be Wigner crystals. 

The reasons for the appearance of Wigner crystal mentioned above also apply to the 
bilayer case. However, because of the additional internal degrees of freedom, bilayer systems 
have a much greater variety of Wigner crystal states. A limited version of this rich variety 
can be seen even without calculations, which we present schematically in figure 1. When 
the layers are far apart, [i.e. regime (i)], the system becomes two single layers, each one 
has its own hexagonal Wigner crystal. To reduce interlayer repulsion, these two crystals 
will be "staggered" as in the usual hexagonal close pack, (i.e. the lattice points of one 
lattice lies directly above the centers of the triangles of the other). As the layer spacing 
D is reduced, one enters regime (ii). In this regime, the staggered hexagonal structure in 
regime (i) cannot survive because the large inter-layer repulsion will push the electrons to fill 
up the "vacancies" in the staggered structure. The simplest staggered crystalline structure 
without vacancies is the staggered square lattice shown in figure 1, where the lattice sites 
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of one lattice sit above the centers of the other. As D decreases further, one enters regime 
(iii). As mentioned before, the bilayer system reduces to the single-layer case. The crystal 
structure is therefore hexagonal. 

The purpose of this paper is to discuss how these crystals transform into one another as 
the energy parameters and the filling factor of the system are varied. We have performed 
a variational calculation which examines a large class of lattice and pseudo-spin structures. 
The resulting phase diagram is surprisingly rich. It contains more crystalline phases than 
the ones mentioned above. It also contains a multicritical point and a critical end-point. A 
key feature of these crystal phases is that their spin structure vary over most parts of the 
phase diagram, even in those regions where the lattice structure remains fixed. In fact, as we 
shall see, the simple-minded staggered square structure mentioned above (which has fixed 
spin structure in the pseudo-spin language) is not the optimal structure in general. Another 
common feature is that all bilayer crystals acquire a net magnetization in the psuedo-spin 
space, except for the case of strictly zero tunneling. The existence of this magnetization 
has important implications on the macroscopic properties of these electron solids, which we 
shall discuss elsewhere. 

The rest of this paper is organized as follows. In Section II, we describe our model and 
our variational wavefunction. The result of our calculation, i.e. the Wigner crystal phase 
diagram, is presented in Section III. In Section IV, we describe in detail our calculational 
scheme. Although our variational scheme is conceptually straightforward, its actually im- 
plementation is not. The reason is that in evaluating the energy of a crystal structure, one 
has to evaluate many lattice sums which are very slowly converging. Moreover, one has to 
repeat the energy calculation thousands of times in the minimization scheme. It is therefore 
necessary to develop tricks to speed up the evaluations of the lattice sums to the point that 
makes the minimization process feasible. When all the tricks are put together, the calcula- 
tion is quite involved. To make the presentation as clear as possible, we first describe the 
computational strategy in detail in Section IV. The details of evaluating the relevant expres- 
sions are given in the Appendices. At the end of Section IV, we also compare our results 
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with the recent work of some others on the same subject. In Section V, we summarize our 
results, and emphasize those important features of bilayer Wigner crystals which we believe 
should be model independent. 



II. THE WIGNER CRYSTALS IN BILAYER QUANTUM HALL SYSTEMS 

The model and the pseudo-spin representation : 

Consider the double-layer quantum well shown in figure 2. We shall adopt the coordinate 
system shown in figure (2a) and (2b). The origin of coordinate system is placed at the center 
of the well. The layers are parallel to the x-y plane. The position of an electron will be 
denoted as x = (r, z), where r = (x,y) denotes a two dimensional vector. We shall consider 
magnetic fields H normal the layers, i.e. along z. The Hamiltonian of a system with N 
electrons is H = H Q + V — B, where H Q is the single particle Hamiltonian, V is the electron- 
electron interaction and B is the (divergent) classical contribution in V (which we shall 
define later). Explicitly, H Q = X^I=i[-^||( r i) + H±(zi)], where H\\ describes the cyclotron 
motion in the xy plane, and H± describes the motion along z, 



m* is the effective mass, and U(z) is a double well potential as shown in figure 2. We shall 
consider the case where the ground state and first excited state of U (z) (denoted as /+ and 
/_ respectively) consists of maxima at D/2 and —D/2. (See figure 2). The separation D 
between the maxima will be referred to as the "separation" between the layers. 
The electron-electron interaction V is given by 

Since the electrons are localized around z = ±D/2, a system will total electron density n 
will contain a "classical" contribution B, which is the Coulomb energy of two infinitely thin 
layers located at ±D/2, each with density n/2, 





2m* dz 2 




(1) 
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|y-y'l 



(3) 



n 



n b (r, z) = - [5(z - D/2) + 5(z + D/2)} 



(4) 



The ground state /+ and the first excited state /_ are symmetric and antisymmetric 
about the center of the well. Their energies will be denoted by e + and e_. For later 
discussions, we define 



A = e_ - e+, 



(5) 



which we shall refer to as the "tunneling gap" . The ground state of H\\ is the lowest Landau 
level with energy \hui C) uj c = eH/m*c, with wavefunctions 



u m (r) = ^MLe-^ 



(6) 



(2n2 m m\i 2 y/ 2 

For large magnetic fields and a sufficiently narrow well, it is sufficient to consider the lowest 
Landau level {u m (x, y)} and the "pseudo-spin 1/2" space spanned by / + and /_. Although 
/ + and /_ are eigenstates of H±, sometimes it is more convenient to use the "localized" 
basis ft and f\ 



/ T = ^(/ + + /-), 



fl — (/+ ~ /-) > 



(7) 



which are states localized in the upper and lower layers respectively, (see figure 2). If we 
define spinors 



/(*) 



c 



\Ci) 



(8) 



a general linear combination of /+ and /_ can be written as (using the fact that / is real) 



/t . c = Ct/t + 



(9) 



which is completely specified by the spinor (. 
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In denning the spinor / in eq.@, we have in fact implicitly chosen a coordinate system 
(xi, x 2 , x 3 ) in the pseudo-spin space so that /j and f [ are eigenstates of a 3 = a ■ x 3 , and / + 
and /_ are eigenstates of o\ — a ■ xi respectively, 



v / 



v 1 / 



u 4, 



(10) 



v- 1 / 



There are no relations between (x!,X2,x 3 ) and the real space directions x, y, and fj_. In 
general, a spinor can be written as 



c 



/ cos | e -^/2\ 



sm~e 



%4>/2 



-ix/2 



(11) 



where 6 and <p are the polar angles of its spin vector S. 



S = C ?C = cos#x 3 + sin^ (cos0 Xi + sin0 x 2 ) . 



(12) 



The spin vectors of / + , /_, fr, and /; are xi, — Xj, x 3 and x 3 respectively. (See fig. 2). In this 
representation, the energy of N non-interacting electrons in the lowest Landau level and in 
the pseudo-spin space is 



H = -\&S U 



N 



Si = S • xi = ^Si,*. 

8=1 



(13) 



where we have ignored the constant N(hu c /2). 



Bilayer Wigner Crystal Variational States : 



In the single-layer case, Wigner crystals can be constructed using coherent states |L0 
A coherent state at R is defined as 



=cxp 



r 
"4£2 



{x + iy){R x - iR y ) K 2 



-4^P <r|R> 



(14) 



A simple variational state for the Wigner crystal can be obtained by antisymmetrizing a 
product of coherent states distributed on a lattice {R«}. To simplify the notation, we write 



(15) 



The Wigner crystal wavefunction is then 



1$™ > = A 



MM)= ^ 



iV 



i=l 
■ N 

n 

i=l 



E(-i) p 
p 

= E(-i) p 



' AT 

i=l 
N 

i=i 



(16) 
(17) 



where [r] = (r 1; r 2 , rjv), .4. is an antisymmetrizer for the electron coordinates {ij}, and 
P denotes a permutation of N objects with signature (— l) p . 

Bilayer coherent states are simply single-layer ones augmented by a pseudo-spin struc- 
ture, 



RjC (x) = (p(z) ■ C) <Mr) x|R, C >, 



(18) 



recalling that x = (r, z), r = (x, y). The analog of the Wigner crystal of the single-layer case 
can be obtained by antisymmetrizing the product of a set of N coherent states, distributed 
on a regular array of N points {Rj, i — 1, ..N} in the x-y plane. Each point Rj is associated 
with a spinor CC^U) describing the wavefunction of this coherent state along r±. As in the 
single- layer case, we write 



(J) = teu.tfROW; 



(19) 



The index % now stands for both the 2D lattice site R, and the spinor C(Ri)> where the bold 
face j denotes the three dimensional coordinate x 3 - of the j-th electron. In this notation, 
Wigner crystal wavefunctions are still of the form eq . ([TBI) . More explicitly, they are 



<£w([x]) = A 



N 



.i=i 



E(-i: 



■ N 

n*«(i) 

.i=i 



(20) 



Since these crystals are completely specified by the two dimensional array {R«} and the 
spinors {£(Rj)}, it is to useful to represent them as "two-dimensional" crystals of "spin 
1/2" particles. In this representation, the hexagonal structure in regime (iii) (i.e. figure 
lc) corresponds to the ferromagnetic state shown in figure 3a, with magnetization along xi, 
representing the symmetric state along r±. The centered square (fig. lb) and staggered 
hexagonal (fig. lc) structures mentioned in Section I correspond to the antiferromagnetic 
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structures shown in figure 3b and 3c. The spins along +x 3 or — x 3 in these figures correspond 
to electrons in the upper or lower layer. 

Ideally, one would like to do a variational calculation taking the set {R«} and the spinors 
{C(Rj)} as variables. Such a parameter space is too large to be practical. To anticipate 
the effects of tunneling on the lattice and the spin structure we have discussed, we restrict 
ourselves to the following configurations: We shall consider systems with an even number 
of electrons, (N even). The array {Ri} consists of two lattices A and B, which are identical 
except shifted relative to one another by a vector c. Spinors on the same lattice are identical, 
but need not be the same on different lattices. In other words, if ai and a 2 are basis vectors 
of A, then 

R= n-ia.i + n 2 a 2 , and C(R) — (a, if R G A, 

R= m iai + m 2 a 2 + c, and C(R) = (b, if R e B. (21) 

where {n^mj} are integers. 

Our variational calculation is performed at fixed electron density n and fixed magnetic 
field. The variational parameters are the spinors (a and (b, the displacement vector c, and 
the basis vectors ai, a 2 . The latter are subject to the constant area constraint |ai xa 2 | = 2/n. 
The energy function to be minimized is 

In the single-layer case, since there is only one energy scale, e 2 /i, the energy per particle is 
of the form 

I = T £ ("»- < 23 > 

where £ is a dimensionless function of v. In the bilayer case, because of the three energy 
scales e 2 y/n, e 2 /D, A, the energy per particle is of the form 

E = j£(D/£,A/(e 2 /£),v) (24) 
10 



The labor of our calculation is to minimize the energy eq.(p2|) within the variational space 
{T = (( A , ( B , c, a 1; a 2 )}. The phase diagram of the system is given by the optimal configu- 
ration T° as a function of the experimental parameters {D/£, A/(e 2 /£), u}. 

In the next section, we shall present the Wigner crystal phase diagram according to our 
calculation. The details of our calculation will be given in Section IV. 

III. PHASE DIAGRAM 

Our variational calculation reveals an unexpected Wigner crystal phase diagram. Be- 
cause of its richness, we shall display it in two different ways. We first show the phase 
diagram in the plane of A/(e 2 /a) and D/a for various filling factors z/, (see figures 4a, 4b, 
and 4c), where 

a 2 = 2/n = 4tt/z/ (25) 

is the area of the unit cell of lattice A. Since the layer spacing D and the tunneling gap 
A are fixed in actual experiments, variations in A/(e 2 /a) and D/a correspond to varying 
the electron density n (say, by varying the gate voltages). Figures 4a, 4b, and 4c are for 
v = 1/3, 1/5, and respectively. Of course, at v = 1/3 and 1/5, the system is a quantum Hall 
fluid. The reason that we still choose to display the Wigner crystal phase diagram at these 
fillings is because it is essentially unchanged at nearby fillings, where the system is no longer 
a quantum Hall fluid. In fact, when plotted in terms of the variables A/(e 2 /a) and D/a, 
the phase boundaries show only relatively small shifts over a large range of filling factors. 
On the other hand, if the phase diagram is displayed in terms of the variable A = A/(e 2 /£) 
and v for different values of D/£, (figures 5a to 5d), the movements of the phase boundaries 
become much more pronounced. Figures 5a to 5d will be useful for experiments where 
electron densities are fixed so that variations in magnetic field causes variations in A and v. 

We have found altogether five stable Wigner crystal states. We shall label them together 
with their region of stability as I, II, III, IV, V. Roughly speaking, the two single-layer 
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regime (i.e. regime (i)) mentioned in Section I is contained in region V. The correlation 
regime (ii) is contained in region II, III, IV. The one-component regime (hi) is contained 
in region I. The optimal spin configuration turns out to be 

Sa(b) = shifts! + (— )cos#x 3 (26) 

The spin structure of crystal I is ferromagnetic, (i.e. 9 = 90°, = Sb = xi). All other 
crystals have mixed ferromagnetic and antiferromagnetic order, with < 9 < 90°. These 
system has a uniform magnetization 2sin#xi and a staggered magnetization 2cos#x 3 . The 
pure antiferromagnetic case 9 = only occurs at zero tunneling. The five classes of stable 
Wigner crystals are : 

Region I : One-component ferromagnet hexagonal crystals : The lattices A and B of these 
crystals are staggered in such a way so that their union is a hexagonal lattice, [a x • a 2 = 
0, a 2 /ai = \f?> ) a\ a/3/2 = 1/n, c = (ai + a 2 )/2]. In the entire region I, Sa = Sb = Xi, i.e. 
all electrons are in the symmetric state in the r± direction. This structure is "rigid" . Both 
its lattice and spin structures are fixed in the entire region I. 

Region II : Mixed ferromagnetic and antiferromagnetic centered rectangular crystals 
: Both sublattices A and B are rectangular lattices, with B sitting at the center of the 
unit cell of A: [ai • a 2 = 0, ai<2 2 = 2/n, c = (ai + a 2 )/2]. Unlike the "rigid" hexagonal 
structure in region I, these structures are "soft" in the sense that both the spin angle 9 and 
the lattice parameter a 2 /ai vary continuously throughout the entire region II. (See figure 
6a). These crystals are separated from the hexagonal ones in region I by a first order line, 
which changes into a second order line at a multicritical point. It is also separated from 
crystal III mentioned below by a second order line, which intersects the I-II first order line 
at a critical end-point. 

Region III : Mixed ferromagnetic and antiferromagnetic centered square crystals : Both 
sublattices A and B are square lattices. They stagger in the same way as those in region 
II. The lattice structure of these crystals are rigid, with [ai ■ a 2 = 0, a\ = a 2 = \j2/n, c = 
(ai + a 2 )/2]. However, the spin structure is soft. The spin angle 9 changes continuously 
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within this region. These crystals are separated from crystal I by a first order line, and 
crystal II by a second order line. 

Region IV : Mixed ferromagnetic and antiferromagnetic centered rhombic crystals : Both 
A and B are rhombic lattices. They stagger in the same way as crystal I and II. 
[ai = a 2 , af sin a = 2/n, c = (ai + a2)/2]. Like crystal II, this structure is soft. The 
spin angle 9, and the lattice angle a between ax and a 2 vary throughout the entire region 
IV. (See figure 6b). These crystals are separated from crystal III by a second order line, 
and crystal V mentioned below by a first order line. 

Region V: Mixed ferromagnetic and antiferromagnetic staggered hexagonal crystals : Both 
A and B are hexagonal lattices. A is sitting above the center of triangle of B. The lattice 
structure is rigid, with [ai = a 2 , a = 60°, a?V3/2 = 2/n,c = (a x + a 2 )/3]. However, its 
spin structure is soft. The angle 9 varies over the entire region V. 

Figures 7 and 8 show the variations of the spin angle along the horizontal and vertical 
line in figure 5a. There is a large change in spin angle (about 20°) across the first order 
line separating the one-component hexagonal and centered square structure. However, the 
change across the centered rhombic and staggered hexagonal first order line is small, (about 
2°). As we shall explain in the next section, the spin angle for all cases is given by : 

f 1 if 7 < A 
sin# = I _ _ (27) 

{ A/7 if 7 > A 

where A = A/(e 2 /£), and 7 is a function of the lattice structure (a 1; a 2 ,c), u, and D/£,. 
For crystals I, III, V where the lattice structure is rigid, (see descriptions I, III, V above), 
7 reduces to 7 , 7 , 7 V which are functions of v and D/£ only. We have plotted these three 
curves in figures 9a to 9d. Readers who are interested in the spin angle 9 in regions I, III, V 



for given D/i can extract their values from eq.(27) using these curves. The determination 
of 9 in regions II and IV is less straightforward, for the lattice structures are soft in these 
regions. One has to first determine the lattice structure by numerical minimization, and 
then evaluate 7 following the prescriptions in Section IV and Appendices C and D. 
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IV. VARIATIONAL CALCULATION 



In this section, we shall describe in detail our variational calculation. Our goal is to 
minimize the energy per particle at fixed densities with respect to the variational Wigner 
crystal wavefunctions of eq.(BO). 



S =< H > /N = [< H > + < V > -B] /N = S + V - B 



(28) 



where < O >=< ^ wc \0\^/ wc > / < ^ wc \^wc >■ To simplify the notation, we shall from 
now on measure lengths in units of magnetic length £, and energy in units of e 2 /£. In these 
units, the tunneling energy is 



So 



1_ 



4 



A(S? + S?), A = A/(e 2 /£) 



(29) 



Separating the background interaction < B > into contributions due to each layer, we have 

(30) 



„ n f d 2 y n 
B=- — + - 



d 2 y 


' 2vr 


71 D 


Vy 2 + n 2 


Go 2 





y 

where a 2 = 2/n = An/v is the unit cell area. Although bilayer Wigner crystals are more 
complex than the single-layer ones, they are identical in form when represented as in eq.(|i~6D 
and eq. (|20|) . The evaluation of < V > can therefore proceed identically as in Maki and 
Zotos |TD[. The result is an expansion in a set of n-body potentials, 



where 



<y >= J2V(ij)+ £ V(ijk) 

i>j i>j>k 



V(ij) 



(31) 



(< ij 


- <ji\)V(\ij > 


>) 


2(1- 


< i\j > 


2 ) 



(32) 



Recalling that i stands for R and C(R-); we have 



1 



1 



d a x / d 3 x' 



2 \l-\S(ij)\ 2 / 

12 J d 3 x 0*(x)0,(x) 



X — X' 



-|0 i (x)0 i (x') - 



X )0i X 



'M2 



IS 



'j i 



(33) 
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The integrals in eq.(|33"D include both direct and exchange terms, which are proportional 
to |0j(x)0j(x')| 2 and 0*(x)0j(x)</>*(x')</>j(x') respectively. As in the single-layer case, the 
exchange terms as well as the n-body term for n < 3 are smaller than the direct terms by 
a factor of exp[— R 2 /2£ 2 ], where R — |Rj — Rj| is the distance between two coherent states. 
This factor is less than 10~ 3 for v < 1/2. Thus, if we focus on Wigner crystals at filling 
factors less than 1/2, the exchange terms in eq.([33"l) and the n-body terms in eq.(|3~Tl) for 
n < 3 can be neglected. V(ij) then takes the Hartree form. 

V(ij) = |d 3 x J dV^-^l^x^x')! 2 , for " < 1/2 (34) 

Next, using 

1_ = r d^g 2vr iq . (r _ r0 K _ rl | 
|x-x'| J (2tt) 2 q v ; 

eq.(|34]) can be written as 

V(ij)= / dg e' q J (qR)W(q) (36) 

J 

W(q)= J dudu' e-^-"'l \f(u) ■ C(R,)| 2 ■ C(R,)| 2 - (37) 

Recall that f\(r±) and fi(r±) are localized around r± = D/2 and —D/2. Typically, these 
functions decay away from the layers within a decay length 1/k, [i.e. f^m (r±) ~ exp(— k \r± — 
(+)D/2|)]. In Appendix A, we show that as long as the products kD and k£ are moderately 
larger than 1, (referred the "moderate" condition), the spinor products in eq . (j37|) can be 
replaced by 

|/t(«) • CI 2 = \Ci\ 2 5(u - D/2) + |C;| 2 % + D/2). (38) 

These "moderate" conditions are certainly feasible in experiments. Within the approxima- 
tion eq. (|38|) , eq. (|36|) reduces to the simple form 

V(ij) = U + (R) + l7_(R)5 8 (i)5 8 (7) = ^w,ay)( R ) ( 39 ) 

C/±(R)= - / dq e-i (1 ± e-o D ) J (qR) (40) 
2 Jo 
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where R = R; - R,, S 3 (i) = S 3 (Ri). 

Note that the ineffectiveness of the exchange in eq.(|33"D and eq.(p?|) does not mean that 
the problem is classical. Quantum mechanical effects are manifested through (a) the smear- 
ing of the (classical) 5-function density to a Gaussian in the lowest Landau level, (thereby 
affecting the interaction of two electrons at distances of the order of £), and (b) layer tun- 
neling, (i.e. the tunneling gap A in eq.(p9|)). [Note that the insignificance of layer exchange 
in eq. (|37l ) under the "moderate" condition does not mean that the tunneling gap in eq.([29|) 
can also be ignored. The reason is that A is of the order of |e|e~ KD , |e| = |e+ + e_|/2. Even 
though e~ nD is small when kD is moderately larger than 1, A can still be comparable with 
with other energies in the system for sufficiently attractive quantum well, which makes e 
sufficiently large.] 

To evaluate V, we separate the contributions from different lattices A and B. Denoting 
S3 and S3 as S and 5" respectively, we have 



v= \ E [Vs,s(R) + Vs>A*)) + \ E VsA* + c ) 

4 R^O 1 R 



7 S£-S 3 
4 I 2 



S\2 



+ 



4 I 2 



fl\2 



7 = 



A= 2 



Ef/+(R) + E^+( R + c ) 

R^O R 

EqR)-EMR + c) 

R^O R 



E^-(R) + E^-(R + c) 

R^O R 



(41) 
(42) 
(43) 

(44) 
(45) 



All the spin dependence in the energy are contained in eq.(^) and eq. (p9|) , which contains 
(S3 , S3) and (Si + Si) respectively. Since both S A and S B are vectors of fixed length, it 
is easy to see that the energy is minimized when S A = S B = 0. Direct plotting shows that 
U-(R) is a positive, monotonic decreasing function of R. This means A is strictly positive. 
7 may be positive or negative depending on whether c is smaller or larger than the shortest 
vector in A. In either case, it is straightforward to show that the optimum spin structure is 
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= -Si = cos9, Sf = S? = sin9 (46) 



Thus the total correlation energy per electron (eq. fl28|) is 



£ = Z+^sm 2 9-±Asm9, £ = rj - - B. (47) 



The optimum spin angle is therefore given by 



f 1 if 7 < A 
sin# = I _ _ (48) 

[ A/7 if 7 > A 



and the corresponding energies are 

(e + |7-|A if 7 <A 

£={ - 2 _ ( 49 ) 

U-A/(4 7 ) if 7 >A 

Thus determination of the correlation energy and the optimum spin structure for a given 
lattice reduces to the evaluation £ and 7. To evaluate these quantities, it is useful to 
introduce the following sums. Let us define the functions 

r°° 2 
F a CR)= / dq e~i J (qR), 
Jo 

roc „ 

F 2 (R) = / dq e- qD e- q J (qR) (50) 
Jo 

and the lattice sums 

Q<= Qi = ^Fi(|R + c|), i = l,2. (51) 

R^O R 

The functions 77 and 7, (hence £ and 7), can now be expressed as 

r l -B=^[Q 1 + Q 1 + Q 2 + Q 2 ]-B (52) 
-l = {Qi-Qi)-{Q2-Q 2 ) (53) 

The evaluation of 77 — £> and 7 (or £ and 7) reduces to the evaluation of the four sums 
{Qi, Q { , i = 1,2}. The reason that we consider the particular combination 77 — B is because 
(as we shall see) each Q term has a divergent (classical) contribution so that their total 
contribution is B. As a result, the combination 77 — B is finite. Likewise, the differences 
(Qi — Qi), 7 = 1,2 are all finite because their divergent contributions cancel each other. 
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A straightforward forward evaluation of these sums proves highly impractical as they 
converge very slowly On the other hand, an essentially exact evaluation is possible if one 
notes the following : (i) The asymptotic forms of Fi(R),i = 1,2 consist only of powers 
in 1/R or l/y/R 2 + D 2 (see Appendix B), (ii) all sums of the form Er^o i?" (2n+1) and 
Y<r]-/(R 2 + D 2 ) n+l l 2 (where n is a positive integer) can be either evaluated analytically 
or transformed into a rapidly convergent series convenient for numerical calculation. (See 
Appendix C). We therefore adopt the following strategy : 

(a) For any desired accuracy, (which we choose to be 10 -9 ), we choose a cutoff A and a set of 
asymptotic functions {F? S (R), i — 1, 2} which are obtained by retaining a few terms of the 
asymptotic expansions of Fi(R) such that {F® S (R)} reduce to {Fi(R)} within the desired 
accuracy for all R > A. (See Appendix B). Next, we rewrite the lattice sums Qi and Q i as 

Qi = QT + QT, Q^QT + QT, (54) 

where 

R^O R 

and 

QT= E [F(R)-Fr(R)]&(A-R) 

R^O 

QT= £ [F,(|R + c|) - Fr(\R + c|)] 0(A - R) (56) 

R 

The superscipt sr means "short range", and Q(x) = 1 or if x > or < 0. It is clear that all 
the short range contributions {QfiQ^} are finite. 

(b) The sums {Q% s , = 1,2} are evaluated analytically by generalizing the method of 
Bonshall and Maradudin |ll[j , or transforming the sum into a very rapidly convergent series. 
These are done in detail in Appendix C and D. As we shall see in Appendix C, all four 
asymptotic sums contain a divergent piece typical of Madelung sums. Separating out this 
divergent piece, they can all be written as 
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where the {Tj,Tj} are the finite parts of {QT iQT}- Using eq.fl52|) and (|53|), £ and 7 in 
eq.(^7|) are now expressed in terms of the finite quantities {Tj, T;, Qf , QT}'- 

v-®= E llQT + QT + Ti + Ti 



7 



[(QT - QT) - (QT ~ Q2) + (Ti - TT) - (T 2 - T 2 )\ (58) 



The explicit expressions of Tj and T, are given in Appendix C. 

(c) The short range sums Qf and Qf can be evaluated directly. This is because the 
asymptotic function F® S (R), is a simple sum of inverse powers of R, and Fi(R) is proportional 
to the zeroth order order Bessel function with imaginary argument, which is available in most 
math libraries. The evaluation of {QfiQT} * s ^ ess straightforward as the integral F 2 (R) is 
not a tabulated special function. Although this function can be evaluated to high accuracy 
by numerical integration, incorporating this integration in the minimization process is too 
time consuming to be practical. However, this bottleneck can be eliminated by the following 
trick. We first evaluate F 2 (R) by numerical integration on a fine mesh within the short range 
region R < A. The values of F 2 (R) at any point not on the mesh can be obtained (to a 
desired accuracy) from the nearby mesh points by using the "cubic spline" interpolation [[T2| . 
This allows us to replace all necessary numerical integrations in our minimization process 
by the spline interpolates, which is enormously faster. 

With the short range contributions given by (c), and the asymptotic contributions given 
in Appendix C and D, our evaluation of the functions £ and 7 is complete. 



The accuracy of the present calculation : The most serious approximation in our 
calculation is that the exchange terms are ignored in eq. (|33|) . As mentioned at the beginning 
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of this section, the exchange terms are smaller than the direct terms by a factor of e _i?2 / 2 ^ 2 , 
which is 1CT 3 around v — 1/2 and drops dramatically at lowering fillings, (1CT 5 around 
v — 1/3 and 1CT 7 around v = 1/5). Once the exchange terms are dropped, the rest of 
our variational calculation is essentially excact. The numerical scheme we have mentioned 
can easily acheive an accuracy of one part in 10 9 , (and can be improved systematically). 
To demonstrate further the accuracy of our calculation, we show in Figure 10(a) and 10(b) 
the energies of vairous crystal states at v — 1/3. The exchange energy is invisible on the 
scale of both figures. One can clearly see on figure 10 that the energy differences between 
different crystals range from 10" 1 (far from the phase boundary) to 10~ 2 (close to the phase 
boundary), except for those between crystal I and II, and between II and III. The latter 
is shown on an expanded scale in figure 10(b), from which one can see that the energy 
difference of these crystals is typically of the order to 10~ 3 . Thus, the energy differences 
between all crystals are much greater than the exchanged energy 10~ 5 . Furthermore, the 
exchange energy appears as a systematic correction to all crystals. In comparing the energies 
between crystals of similar crystal structures, they are to a large extent cancelled out so 
that the actual exchange contributions to the energy difference are at least an order of 
magnitute smaller than 10~ 5 . The only place where our calucaltion may begin to produce a 
few percent error bar is close to v — 1/2. At lower fillings, all numerical evidence indicate 
that our evaluations of crystal energies are accurate to the order of e^l v ', (which is the 
ratio between exchange and direct terms), and that the energy difference has an even higher 
accuracy (10~ 1 e~ 7r / Jy ) because of the aforementioned systematic correction. 



V. CONCLUDING REMARKS : 

We have seen that Wigner crystals in bilayer quantum Hall systems come in different 
"magnetic" and structural varieties. Although our conclusions are based on a variational 
calculation, there are two characteristics of these crystals which are direct consequences 
of the competition between layer tunneling and interlayer Coulomb repulsion, and should 
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therefore exist in real systems. They are the coexistence of antiferromagnetic and ferromag- 
netic order, and the coupling between lattice and spin structure. The antiferromagnetism in 
pseudo-spin space is to minimize the interlayer Coulomb repulsion, whereas the ferromag- 
netism is to take advantage of the tunneling energy. The coupling between lattice and spin 
configurations simply reflects the competition of these two energies. In many ways, Wigner 
crystals in bilayer systems are like 3 He solids, which also have different phases with different 
lattice and magnetic structures. The magnetic structures in 3 He solids are known to give 
rise to many remarkable nuclear magnetic resonances. It will not be surprising if bilayer 
Wigner crystals also have many interesting "magnetic" properties. 
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Appendix A. Conditions for eq.(|38|) : 

The integrand in eq.([37|) contains terms like f-\(u) 2 fi(u') 2 , f^(u)fi(u)fi(u') 2 , 
fl(u) 2 f^(u') 2 , etc. We shall discuss the contributions of the first two terms, from which 
the contributions of the rest can be similarly inferred. To simplify these expressions, we 
shift the origins of u and v! so that the maxima of /f and f\ are at the origin, 
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The contribution of fa (u) 2 fi(u') 2 and fa(u)fi(u)fi(u') 2 to eq.(|37|) can then be written as 

h{q) = JdsJ ds'e^ D+s - s '\g{sfg{s') 2 (60) 

I 2 {q) =JdsJ ds'e^ D+s -^g(s)g(s + D)g(s') 2 (61) 

Since g(s) ~ e~ K ' s ', the exponential factor in eq.(|5U|) and fl5"T| ) can be replaced by e -i( D + s - s ') 
if kD is moderately larger than 1. On the other hand, the Gaussian factor in eq. (|36|) limits 
q to the range q < l/£. Hence, if k£ is moderately larger than one so that e~ qs decays much 
slower that e~ KS , the functions g{s) 2 and g(s') 2 in eq.(|60D will act like ^-functions. The same 
approximations applies to 12- However, because of the product, I2 is smaller than I\ by a 
factor of e~ KD . It can therefore be ignored if kD is moderately larger than 1. Applying the 
analysis to other products of /'s, we have eq. (|38|) . 



Appendix B : Asymptotic expansion of F\(R) and -^(-R) 



Note that 



F 1 (R) = / dqe-i J (qR) = ^-1^/8). 
Jo 2 



(62) 



where Il{x) = e x I Q (x), and I D is the Bessel function with imaginary argument. For R ^> 1, 
Fi has an asymptotic expansion 



~ 1 1 00 \(2n - lVH 2 1 
F 1 (i?)-F 1 -(R) = E i 7V 2 4 = E- ' 



n=0 n - 11 n=0 



(63) 



n\ R 2n+l 

For R > A = 35, (i.e. A = 35£ in the original unit), only three terms in the asymptotic 
series 



1+V 2 + -V 4 



F? S (R) 

L ^ 

is enough to reproduce -Fi(-R) to the accuracy of 10~ 9 



1 

R 



(64) 
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In the case of F 2 (R), we note that it has the asymptotic expansion 



f 2 (r)^f™(r) = y: 



n=0 



-1 



n\ dD 2r > 



VR 2 + D 2 ' 



(65) 



£( . ir M 

71=0 ' 4 - 



! (i? 2 + £)2)n+l/2 



2n 



VWTd 2 



(66) 



P n (x) are the Legendre polynomials. As in the case of Fi, for R > A = 35, only three terms 
in the above series 



F?(R) 



1 d 4 



1 



V^ 2 + D 2 



(67) 



d£> 2 2dL> 2 _ 

is enough to reproduce Fi(R) to the accuracy of 1CT 9 . P n (x) are the Legendre polynomials. 
The functions F^ S (R) and F£ S (R) are the asymptotic functions we have used in our numerical 
calculation. 
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Appendix C : Explicit determination of the asymptotic lattice sums {Q°; s ,Q^ s ,i = 1,2}: 
Consider the following generating functions : 

1 



He) =£■ 



R 



R+c| 



(68) 



^ (\R + c\ 2 + D 2 y/ 2 

where Re A. Using the asymptotic functions {F® s (R) , F£ s (R)} defined in Appendix B, it 
is easy to see that 



Qf = Lim b ^o 



1 + V£ + -V 4 



1, 

2 



h(b) 



m 

Q 

q7 



as _ 
2 " 



V" + -V 4 
v c 1 2 



h(c) 



d 2 




1 d 4 " 




+ 




dD 2 


2dL> 4 


d 2 




1 d 4 : 




+ 




dD 2 


2dL> 4 



9(0, D)-- 



g(c,D) 



(69) 



We shall first give the expressions of h and g. Their derivations are given in Appendix D. 
These expressions are 

2tt 



He) 

g{c,D) 



Ga 2 
2tt 



Ga 2 



G-»0 



+ - [-2 + L(c/a)} [ L (I R + c l/ a ) + L(R/ a)cos(2n\R x c|/a 2 )] (70) 



2nD 1 



ar a 



R^O 

e -2wRD/a 2 

£ rwi \ cos ( 27r l R x c l/ fl2 ' 

K ^ [R/a) 



(71) 



where a 2 is the unit cell area of lattice A, (a 2 = 2/n = 4ml 2 jv), and 



l(x) = - (i - <KV*x)) , <P(v) = 4= f du e ~" 

X v 7 \ 71 JO 



(72) 

y /i j\j 

The function is the error function. 

It may seem that only eq. (|7T|) is necessary because h(c) is a limiting case of g(c, D). The 
latter is true but impractical. The reason is that the number of terms needed to include in 
the sum in eq.(|7T|) to achieve a specified accuracy grows as 1/D as D — > 0. The series in 
eq.flTTD is therefore useless in the small D limit. We are, however, lucky for two reasons. 
First of all, the typical value of D in real experiments is of order unity. The series expansion 
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eq.(|71~D is highly convergent for these D values. As for h(c) [which corresponds to the special 
case g(c,D = 0)], an analytic expression [eq.([T(])] can be obtained by a straightforward 
generalization of the method of Bonshall-Maradudin (BM) ||11|| . (See appendix D). This 
method produces a super convergent series for h, (a cutoff of just 4 lattice constants in the 
sum is usually sufficient to produce an accuracy of 10~ 12 ). 

Next, we note that the gradients in eq.(^) can be conveniently evaluated using the 
following identities: 

V 2 L(r) =L 2 (r) = \L(r) + pi(nr 2 )e~ nr2 ] /r 2 , p 1 {x) =Ax + 2 

V 4 L(r) = L 4 (r) = [9L(r) + p 2 (nr 2 )e~ 7rr2 ] /r 4 , p 2 (x) = 16x 3 - 8s 2 + I2x + 18 (73) 

With these identities, using the definitions in eq.(|57|), and eq.(|69|) it is straightforward to 
work out the finite parts {Tj,Tj,z = 1,2} of the asymptotic sums {Q™, , i = 1,2}. For 
Qi s , we have 

2tt 



Lim b ^ (h(b) - 1/6) 

72 



G^O 
16 



Go? 

Um.^ Vi{h{h)-l/b)=T, 
^Lim b ^oV 4 (h(b) - 1/6) = T u 
T\ = Ti a + Ti b + Ti c 



+ Ti a 



(74) 



where 



Tl a =- 

a 



Tlc "2^ 



-4 + 2 L(R/a) 

R^O 

+ E (HR/a) - Air 2 (R/a) 2 L(R/a 



R^O 

64vr 2 



+ (U(R/a) + 16ir 4 (R/a) 4 L(R/c 

R^O 



(75) 



Similarly, for , we have 
h{c) 



Ga 2 



+ T la , V 2 c h(c)=T 



ib- 



G->0 

Ti = T Xa + + Tic 



\v A Mc) = Ti c 



(76) 



25 



where 



T - 1 

J- la— - 

a 



T - 1 

J- lb— — ? 



L(c/a) + ( L (I R + c l/ a ) + L{R/a) cos(2tt|R x c| 



R^O 



L 2 (c/a) ( L 2(|R + c|/a) - An 2 (R/a) 2 L(R/a) cos(2tt|R x c|/a 2 )) 



R^O 



lc- 



2a 5 



L 4 (c/a) + ]T (L 4 (|R + c|/a) + 16n A (R/a) A L(R/a) cos(2tt|R x c|/a< 



R^O 



Next, for Qf, we have 



[g(0,D)-l/D] 
d 2 



2tt 



dD 2 
1 d 4 



(g(0,D)-l/D) = T 2b , 
(g(0,D)-l/D)=T 2c 



2nD ) + T 



2dL> 4 
T 2 = T 2a + T 2 6 + T 2c 



where 



-2a = 



_j _ e -2nRD/a 2 

(^) + rVo R / a 



T 2b = 



2! 



T 2c - 



1 

2^5 



(D/o) 
4! 



+ ]T 47r 2 ( J R/a)e- 2wDi? / a2 



R^O 



+ l^\R/afe-^ DR ' a2 



R^O 



Finally, for Q 2 , we have 



g(c,D) 
d 2 



2tt 



dD 2 
1 d 4 



»(c,o) = r a , 



h~ + J 2a, 

a 2 / 



2dD 4 ' 
T 2 = T 2a + T 2 b + T 2c 



where 
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J e -2irRD/a 2 

a r^o W a ) 



T 26 = V 4ir 2 (R/a)e~ 2nDR/a2 cos(2vr|R x c|/a 2 ), 

— ' 

" R=£0 

T 2c = £ 16vr 4 ( J R/a) 3 e- 2 ^ /a2 cos(27r|R x c|/a 2 ). (81) 

^ a R^O 

A final note — to obtain an overall accuracy in the correlation energy to eight significant 
digits, suitable cutoffs (Ah and A g ) for the lattice sums in eq.([70D and (|7lD are 5a and 
A g « 20a 2 /7rD (so that e" 27rDA ^ a2 ~ KT 11 in eq.flTI])). 



Appendix D: Derivation of h(c) and g(c,D). 

The series for /i(c) in eq.([70D can be derived by a straightforward generalization of the BM 



method [11|. We first write 

Kc) = E ^ (jf° + jT) <tf e^ 2 ' R + c ' 2 ee (0 + («) (82) 

where w = y/n/a. After rescaling the integration variable, the second term (ii) can be 
written as 

(«) = i$>(l R + c l/°). (83) 

a R 

where £(-R) is defined in eq . (|72|) . To evaluate the first term (i), we first convert the real 
space lattice sum into a sum in the reciprocal space, (using the standard relation Sr^(R) = 
X)g^(G)/ci 2 , where u is the Fourier transform of u), 



Changing f3 — > w/j3, we have 

2 



[t = - 
a 



£ d/3 e -(G/3/2 W )^+iG.c = J2l G (85) 
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The G = term is connected to the 1/R divergence of the sum and has to be treated 
separately: 



^g=o — -LiniG^o 
a 



OO rl 



Jo 



d(3 e 



-(Gf3/2w) 2 +iG-c 



2lT 



Go 2 



2 
a 



(86) 



For the G ^ terms, we note that reciprocal lattice vectors G and real space lattice vectors 
R are related as G = 27rz x R/a 2 . We can then write 



Eg ^ i g = I £ 2 f°°d(3 e -^ 2 P 2 /«^\^\/«* 

a R^O Jl 



L(R/a) cos(2vr|R x c\/a 2 



R^O 



Substituting eqs.(|33|), fl36|), and (|S7|) into eq. (|S2|) we have eq. (|70| 



(87) 



The derivation of eq.flTTD for g(c, D) is as follows. We write g(c, D) as 



Using the relation 



(88) 



we find 



$>^ R = £(2^(q-G)A 

R G 



2tt 



iGc-GD 



(89) 



(90) 



which is eq. (|7T|) if we treat the G = term separately, and and write G = 27rz x R/c 
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Figure Captions 

Figure 1 : (la) Staggered hexagonal structure in regime (i); (lb) Centered square structure 
in regime (ii); (lc) Single component hexagonal structure in regime (iii). The electrons are 
in the symmetric state of the double quantum well, which are represented schematically by 
a "peanut" shape. 

Figure 2 : A simple model of the double layer quantum well and the pseudo-spin represen- 
tation for the lowest doublet. The coordinate system and the quantum well potential are 
shown in fig. (2a) and (2b). The ground state /+ and the first excited state /_ are shown 
in fig. (2c) and (2d). They are represented by spinors (\) and (^J respectively, which have 
spin vectors along +X4 and — X4. Fig.(2e) and (2f) represent the sum and difference of these 
states, which are localized on the right and left well resp. These states are denoted as j\ 
and fi, and correspond to the spinors (q) and (5). The corresponding spin vectors are x 3 
and — x 3 . 

Figure 3: The pseudo-spin representation of the crystals shown in figure 1. The one- 
component hexagonal crystal in Figure lc is represented by the ferromagnet hexagonal 
crystal in figure 3a, with magnetization along Xi. The centered square structure in figure 
lb is represented by an antiferromagnet crystal in figure 3b, with sublattice magnetization 
along +X3 and — x 3 , which are represented by solid and empty circles respectively. This 
structure only exist at zero tunneling A = 0. When A ^ 0, both up and down spins will 
tilt toward xi as shown in figure (3b'). An arrow with a solid or (empty) circle attached 
to the base denotes a spin with positive (negative) X3 component. The staggered hexagonal 
structure in figure la is represented by the antiferromagnetic structure in figure (3c). Like 
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(3b), (3c) only exists at zero tunneling. Nonzero tunneling will make the spins point along 
Xx as shown in (3c'). For later use, the structures in (3a), (3b'), and (3c') are denoted as I, 
II, and III respectively. (See figure 4a to 4c). 

Figure 4 : Figure (4a), (4b), and (4c) show the Wigner crystal phase diagram in the the plane 
of A/(e 2 /a) and D/a at filling factors v = 1/3, 1/5 and 0, where a 2 is the area of the unit 
cell in lattice A. Although the system is a quantum Hall fluid at these fillings, the Wigner 
crystal phase diagram remains essentially the same at nearby fillings, where the system is 
no longer a quantum Hall fluid. The solid and dashed lines are first and second order lines 
respectively. The label (F) means ferromagnetic order, (M) means mixed ferromagnetic 
and antiferromagnetic order. Both spin and lattice structure undergo discontinuous change 
across the first order line. The lattice structure in I, III, and V are "rigid", in the sense 
that they are unchanged in the entire I, III, and V region. In region II and IV, the lattice 
structure (i.e. ai,a 2 , and c) varies with A/(e 2 /a) and D/a. The spin structure is only 
"rigid" in region I where it points along x x . In all other regions, the spin angle varies with 
A/(e 2 /a) and D/a. See Section III for the determination of the spin angle 9. 

Figure 5 : Figure 5a to 5d show the Wigner crystal phase diagram as a function of A = 
A/(e 2 /£) and v for D/£ = 3, 2, 1, 1/2. When plotted in these variables, the phase boundaries 
show more changes when compared with those in figure 4a to 4c. In figure 5a, both phase 
II and the multicritical point are not shown for they appear at very small A and v. They 
are shown in figure 5b, 5c, and 5d. The spin angle 9 along the vertical and the horizontal 
line are shown in Figure 7 and 8. 

Figure 6 : Figure 6a and 6b show the centered rectangular and centered rhombic structure, 
(i.e. crystal II and IV) in figure 4a to 4c. The meaning of the arrows is identical to that 
in the caption of figure 3. The crystal (fig. 6a) reduces to that in fig. 3a when a 2 /a 1 = \/3, 
and 9 = 90°. 
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Figure 7 : The spin angle 9 along the vertical line in figure 5a. There is a close to 20° 
discontinuity across the I-III first order line, where as a very small discontinuity (almost 
invisible on this scale) across the IV-V first order line. 

Figure 8 : The spin angle 9 along the horizontal line in figure 5a. The discontinuities of the 
spin angle across the I-III and IV-V first order lines are similar to that in figure 7. 

Figure 9 : Figure 9a to 9d show 7 I ,7 III ,7 V as a function of v for D/£=3, 2, 1, 1/2. Using 
eq . (f48|) , these curves allow one to determine the spin angle 9 for the crystals I, III, and V 
in figure 5a to 5d. 

Figure 10 : Figure 10(a) shows the energies of the crystal states I to V at v = 1/3 and 
A/(e 2 /a) = 0.1 as a function of layer separation D/a. Figure 10(b) shows the transition 
region I-II and II-III in figure 10(a) on an expanded scale. The energy difference depicted 
is typically of the order of 10 -3 . 
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